#### Simulate BR term structures

# Simulate rates out tt=500 years
N=10e3
tt=500

#### Sigma's posterior quantiles, from Michael Bauer:
#     sigma_u:
#  1st Qu.  Median   3rd Qu.  
#   0.1946  0.2039   0.2140  
set.seed(1)
# Use inverse gamma distribution (the BR prior distribution), calibrated to match
# the provided posterior quantiles
adj.factor = (1.003*0.2047/0.10)^2 # this adjusts the inverse gamma to match the BR posterior
sig = sqrt(1/rgamma(N, shape=100/2, scale=0.04*(100+2)/2)*adj.factor)

quantile(sig, probs = c(0.25, 0.5, 0.75)) # should come close to 0.1946, 0.2039, 0.2140

# rescale so it's not in percentage points anymore
sig = sig/100 

# Generate innovations/shocks
set.seed(1)
head(rnorm(tt*N, 0, sig)) # check values

eps = matrix(NA, nrow=tt, ncol=N)
for (t in 1:tt) eps[t,] = rnorm(N, 0, sig)
eps[1,] = 0 # first period has no shock to fix starting point

# Accumuated shocks
eps.cum = apply(eps, 2, cumsum)

# Set up array to store results
ce.store = array(NA, dim=c(tt, 7),
                 dimnames=list(1:tt,
                               c('1%','1.5%','2%','2.5%','3%','4%','5%')))
r0s = c(0.01, 0.015, 0.02, 0.025, 0.03, 0.04, 0.05)

for (i in 1:length(r0s)) { # for each starting rate...
    # ...simulate real rate, constrained st avg rate>0 ----
    r0 = r0s[i] # starting rate
    
    r = r0 + eps.cum
    
    # We're imposing the constraint that average rates cannot be negative.
    # This is the same as saying the discount factor cannot exceed 1, DF>=1, 
    # so that we get r^ce_t = log(DF)/-t >= 0
    # Do this by calculating the DF, then flooring it:
    DF = apply(r, 2, function(x) exp(-cumsum(x)))
    
    DF = pmin(DF, 1)
    
    # Check:
    # Check forward rates:
    r.floor = matrix(NA, nrow=nrow(DF), ncol=ncol(DF))
    r.floor[1,] = r[1,]
    r.floor[2:tt,] = - (log(DF[2:tt,]) - log(DF[1:(tt-1),]))
    
    # Check time average rates: raw r versus floored
    r.avg = apply(r, 2, function(x) cumsum(x)/seq_along(x))
    r.floor.avg = apply(r.floor, 2, function(x) cumsum(x)/seq_along(x))
    
    range(r); range(r.floor) # note: forward rates can be negative, even after flooring...
    range(r.avg); range(r.floor.avg) # ...but time average rates are >=0 after flooring
    
    E.DF = apply(DF, 1, mean)
    
    ce = log(E.DF)/-c(1:tt)
    
    ce.store[,i] = ce
}
rm(list=c('DF','eps','eps.cum','r.avg','r.floor','r.floor.avg'))

ce.BR = ce.store
save('ce.BR', file=data.dir%&%'/BR_term_structures.RData')
